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Abstract Some novel strategies have recently been proposed for single hidden 
layer neural network training that set randomly the weights from input to hidden 
layer, while weights from hidden to output layer are analytically determined by 
pseudoinversion. These techniques are gaining popularity in spite of their known 
numerical issues when singular and/or almost singular matrices are involved. In 
this paper we discuss a critical use of Singular Value Analysis for identification of 
these drawbacks and we propose an original use of regularisation to determine the 
output weights, based on the concept of critical hidden layer size. This approach 
also allows to limit the training computational effort. Besides, we introduce a 
novel technique which relies an effective determination of input weights to the 
hidden layer dimension. This approach is tested for both regression and classifi¬ 
cation tasks, resulting in a significant performance improvement with respect to 
alternative methods. 

Keywords pseudoinverse matrix • weights setting • regularisation • supervised 
learning 


1 Introduction 

The training of one of the most common neural architecture, the single hidden layer 
feedforward neural network (SLFN) was mainly accomplished in past decades by 
methods based on gradient desce nt, and among them th e large family of tech¬ 
niques based on backpropagation (|Rumelhart et ail Il986ll . The start-up of these 
techniques assigns random values to the weights connecting input, hidden and 
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output nodes that are then iteratively modihed according to the error gradient 
steepest descent direction. Some common drawbacks with gradient descent-based 
learning are anyway the high computational cost because of slow convergence and 
the relev ant risk of convergin g to poor local minima on the landscape of the error 
function ( LeCun et al 

The idea of using the simple and efficient training algorit hms of radial basis 
funct ion neural networks, based on matricial pseudoinve rsion (IPoggio and GirosI 
Il990ll . also for SLFN learning was initiall y suggested in (l^ancelhera 2001); some 
appealing technique s were than developed ( Neuven et all 1201(11 : Kohno et al.L 2011)1: 
A jorloo et fflZ.I . [20(^ and among them the extreme learning machine ELM, (jlduang et al\ . 


20061) which has been successfully applied to a number of real-world applica 


a pplica - 

I I2OIOI) , 


tions ( Sun et ^ . bOOSHWang et al.l . l200SlMalathi et al.l . l2010l : [Minhas et al. 

showing a good generalization performance with an extremely fast learning speed. 

ELM main result states that SLFNs with randomly chosen input weights and 
hidden layer biases can learn distinct observations with a desired precision, pro¬ 
vided that activation functions in the hidden layer are inhnitely differentiable. 

After input weights and hidden layer biases have been randomly set, output 
weights are directly evaluated by Moore-Penrose generalised inverse (or pseudoin¬ 
verse) of the hidden layer output matrix: these two steps conclude one training 
phase and weights are no more modihed, so that their determination is no more 
iterative in the sense of back-propagation based techniques. 

Besides, all pseudoinversion based methods are multi-start, i.e. the above pro¬ 
cedure is repeated meny times in order to hnd a good minimum of the error surface. 
Each training procedure so implies many random settings of input weights and as 
many evaluations of output weights through pseudoinversion. 

However, such techniques seem to require more hidden units than typical values 
from back propagation training to achieve comparable accuracy, as discussed in Yu 
and Deng ( Yu and Pend . [ 2 OI 2 I) . Moreover, pseudoinversion, commonly evaluated 
by Singular Value Decomposition (SVD), is a powerful method but some caution 
is required, since its numerical instability is a well known issue when singular and 
almost singular matrices are involved. 

One aim of this paper is the analysis of these instability iss ues; a preliminary 


assessment of the context and our initial results are discussed in ( Cancelliere et al. 


I 2 OI 2 I) . Here we present further advances on the theoretical framework and we pro¬ 
pose a novel approach to carry out a more efficient learning, showing how singular 
values of SVD can be used to detect the occurrence of numerical instability. 

Besides, we prove the existence of a critical hidden layer dimension that allows 
a careful tuning of the regularisation parameter and the use of regularisation to 
replace unstable, ill-posed problems with well-posed ones. We also propose an 
original method to set input weights that links their size to the hidden layer 
dimension and we show its effectiveness. 


In section[2]we introduce the notation used for describing SLFN architectures, 
and the main ideas concerning input weights setting and output weights evaluation 
by pseudoinversion. In section [3] we discuss the problem of ill-posedness, and the 
basic regularisation concepts. 

In section|3]our framework is tested on some applications selected from the UCI 
database. A substantial improvement in performance with respect to unregularised 
state-of-the-art techniques is shown. 
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2 Input and output weights determination 

Fig. [1] shows a standard SLFN with P input neurons, M hidden neurons with 
non-linear activation functions 0, and Q output neurons with linear activation 
functions. 



Fig. 1 A Single Layer Feedforward Neural Network 


If we have a training set made by N distinct training samples of (input, output) 
pairs (xj , t j ), where Xj £ and tj £ , the training aims at obtaining the matrix 

of desired outputs T £ when the matrix of all input instances X £ R'^^^ is 

presented as input. 

We emphasise that, in the state of the art pseudoinverse approach, input 
weights Cij (and hidden neurons biases) are randomly sampled from a uniform 
distribution in a fixed interval and no longer modified. Therefore this step gives 
the actual input weights values. 

After having fixed the input weights matrix C, the use of linear output units 
allows to determine output weights Wij as the solution of the linear system 


HW = T, 


( 1 ) 


where H £ R 


NxM 


is the hidden layer output matrix of the neural network. 


H = ${XC). It is important to underline that because usually the number of 
hidden nodes is much lower than the number of distinct training samples, i.e. 

M << N, H is a rectangular matrix. 

_ T he least-squar es solution W of the linear system ([1]), as shown e.g. in ( Penrose and To^ . 

ll956l:lBishol3.l20n(ili . is IF = H~^T, where //"*“ is the Moore-Penrose generalised in¬ 
verse (or pseudoinverse) of matrix H. It can be computed in a comp utationally sim- _ 

pie an d accurate way by using singular value decomposition (SVD) fColub and Van Loanl 

We know that every matrix H £ can be decomposed as follows: 


H = UEV^ 


( 2 ) 


where U £ and V £ are orthogonal matrices and X £ is 

a rectangular diagonal matrix whose elements = Cj, called singular values. 
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are nonnegative (usnally the singnlar values are listed in descending order, i.e. 
(Ti > (T 2 > • • • > o-p > 0, p = min {A^, M}, so that E is nniquely determined). 

The psendoinverse matrix has the form 

H+ = VE+U'^, (3) 


where E^ is obtained from E by taking the reciprocal o f each non-zero element 
CTi, and transposing the resulting matrix ( Rao and Mitral Il97lll . The presence of 
very small element oi is therefore a potential drawback of this method. 

For compntational reasons, the elemen ts Oj eqnal to zero or smal ler than a 
predehned threshold are replaced by zeros ( Golnb and Van Loanl [TflQfih . 


3 Singular value decomposition of regularised problems 

To turn an original ill-posed problem into a well-posed one, i.e. roughly speaking 
into a problem insensiti ve to small changes in train ing conditions, regularisation 
methods are often used (Badeva and Morozov. 19911). and Tikhonov reEularisation 

rin[[^ -- 


is one of the most common ( Tikhonov and Arseni m Il977t iTikhono^ Il963ll . 


The error functional to be minimised is characterized by a penalty term Eji 
that depends from the so called Tikhonov matrix F: 


e = Ed + Er = \\hw - T\\l + \\rw\\l, 


(4) 


This matrix can for instance derive from the choice of using highpass operators 
(e.g. a difference operator or a weighted Fourier operator) to enforce smoothness. 
The regularised solution, that we denote by IT , is now given by: 


w = {h'^h + r'^ry^H'^T. 


(5) 


The penalty term improves on stability, making the problem less sensitive to 
initial condit ions, and contains model c omplexity avoiding overhtting, as largely 
discussed in ( Gallinari and CibasL Il999|l . 

To give preference to solutions W with smaller norm ( Bishool . l2flO(T|l a frequent 
choice is F = %/A7, so eqs. © and m can be rewritten as 


E = Ed + Er = \\HW - TWl + X\\W\\l 


w = (h'^H + Xiy^H'^T. 


( 6 ) 

(7) 


The role of the control parameter A is to trade off between the two error terms 
Ed and Er. If A = 0, eq. © reduces to the unregularised least-squares solution, 
provided that {H'^F[)~^ exists. 


_ The regularised solution © can also be expressed (see e.g. (|Fuhrv and Reich^ 

l2012lB as: 

W = VDU'^T ( 8 ) 

where V,U are from the singular value decomposition of H (eq.©) and F is a 
rectangular diagonal matrix with elements 


D, = 


+ A 


(9) 
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Table 1 The UCI datasets and their characteristics 


Dataset 

Type 

N. Instances 

N. Attributes 

N. Classes 

Abalone 

Regression 

4177 

8 

- 

Cpu 

Regression 

209 

6 

- 

Delta Aileron 

Regression 

7129 

5 

- 

Iris 

Classification 

150 

4 

3 

Diabetes 

Classification 

768 

8 

2 

Landsat 

Classification 

4435 

36 

7 


obtained using the singular values of H. 

It is evident that, when unregularised pseudoinversion is used, the presence 
of very small singular values can easily causes numerical instability in H'^\ on 
the contrary, regularisation has a dramatic impact because, even in presence of 
very small values ai of the original unregularised problem, a careful choice of the 
parameter A allows to tune singular values Di of the regularised matrix, preventing 
them from divergence. 

It is clear at this point that a suitable value for the parameter A has to derive 
from a compromise between the necessity to have it sufficiently large to control 
the approaching to zero of a in eq. © while avoiding predominance of penalty 
term in eq.(|6]). Its tuning is therefore crucial to simultaneously control numerical 
instability and overfitting. 

In the next section we propose a strategy to obtain this result showing that we 
achieve better performance and more stable solutions. 


4 Experiments and Results 


Some numeric al instability issues hav e already been evidenced in our previous 
investigations I Cancelliere et al ]. l2ni2|l : we provided suggestions on possible mit¬ 
igation techniques like selection of a convenient activation function and normali¬ 
sation of the input weights. Hereafter we show that adding regularisation to the 
implementation prescriptions already analysed provides a convenient and effective 
approach to deal with such problem. 

The use of sigmoidal activation functions has recently been subject of debate 
because they see m to be more easily driven towards saturation due to their non¬ 
zero mean value ( Bengio and Glorotl I 2 OIOII . while hyperbolic tangent seems less 
sensitive to this problem. 

We select therefore both activation functions for a test aiming at comparing 
their performance in a context where we also compare our proposed regularised 
approach and the unregularised one. Four different experimental settings will be 
analysed, namely HypT-reg, Sigm-reg, HypT-unreg and Sigm-u nreg. _ 

T o further mitigate saturation issues, in our previous work ( Cancelliere et all 
l2012ll we selected input weights according to a uniform random distribution in 
the range {—1/y/M, I/v^M), where M is the number of hidden nodes. This links 
the size of input weights, and therefore of hidden neurons inputs, to the network 
architecture, thus forcing the use of the almost linear central part of the hyperbolic 
and sigmoidal activation functions when exploring the performance as a function 
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of an increasing number of nodes. Such prescriptions are retained in the current 
work. 

We emphasise that so doing input weights are automatically chosen ” small”, 
because the size of interval (—1/%/M, 1/%/M) decreases when the number of hidden 
neurons M increases: for instance with 10 hidden neurons, weights values are 
roughly selected in the range (—1/3,1/3), with 100 hidden neurons in the range 
(—1/10,1/10), and so on. 

Because of the wide use among researchers belong i ng to ELM-community (see 
for instance ( Helmv and Rasheedl . l200^ iHuang et al . 1200(11 : ISun et all . 12001^ our 
performance is also compared with that from unregularised pseudoinversion, input 
weights selected according to a random uniform distribution in the interval (—1,1) 
and sigmoidal activation functions (hereafter, ELM). 

The numerical experiment compares these frameworks applying them to six 
benchmark datasets from the UCI repository ( Bache and LichmanL l^Ol^ . listed 
in Table [TJ 

For each proposed method the number of hidden nodes is gradually increased 
by unity steps, and, for each selected size of SLFN, average RMSE or average mis- 
classification rate on the validation set are computed over a set of 100 simulation 
trials, i.e. over 100 different initial choices of input weights. All simulations are 
carried out in Matlab 7.3 environment. 

Figure [2] gives an insight on the performance trend (resp. average RMSE for 
regression or average misclassification rate for classification tasks) as a function 
of hidden space dimensionality for the cases HypT-reg, HypT-unreg and ELM. In 
Figure [3] performance trends are shown for the cases SigmT-reg, Sigm-unreg and, 
for the sake of comparison, again ELM; because of their similarity with the HypT 
cases, plots are shown for only two datasets (i.e. Abalone and Iris). 

It is interesting to note that when unregularised techniques are used, all datasets 
except Landsat show a fast error growth; besides, the curves have different char¬ 
acteristics for HypT-unreg and Sigm-unreg on one side and ELM on the other, 
showing error peaks in the formers while monotonically increasing error values 
are obtained in the latter. We conjecture the presence of two distinct phenomena: 
numerical instability and overfitting. 

In order to address the former issue, i.e. numerical instability, we evaluate 
for each dataset the ratio between the minimum singular value of hidden out¬ 
put matrix H and the Matlab default threshold (below which singular values are 
considered too small and therefore treated as zero). 

We checked that for each dataset processed with HypT-unreg or Sigm-unreg 
methods there is a critical hidden layer size above which the ratio becomes smaller 
than one; its trend is plotted in logarithmic units (red line, right scale), in Figure 
[2] for HypT-unreg case. 

When approaching critical size, inversion of singular values causes a wrong 
evaluation of and therefore a significant growth in the error; when the crit¬ 
ical dimension is reached, singular values under threshold are automatically re¬ 
moved, thus allowing the subsequent decreas e of error. The same tr end was de¬ 
tected analysing the astronomical dataset in ( Cancelliere et ai] . l2012ll . 

This decrease is anyway not sufficient to reach optimal error values because 
of overfitting, which is known to arise when a large amount of hidden neurons is 
available to reproduce almost exactly the training data. 
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Number of hidden neurons 


Number of hidden neurons 




Fig. 2 Performance comparison in the HypT (-reg and -unreg) cases: Abalone (Top Left); Cpu 
(Top Right); Delta Ailerons (Mid Left); Iris (Mid Right); Diabetes (Bottom Left); Landsat 
(Bottom Right). 


An even more severe overfitting affects ELM in fact in this case test error is 
a monotonically increasing function of the number of hidden neurons. A possible 
explanation is that the setting of input weights in the interval (-1, 1) may al¬ 
low ‘specialisation’ of some hidden neurons on particular training instances, thus 
creating a sort of network partition, carried out thanks to saturation. On the con¬ 
trary, when weights are randomly selected in the interval (—1/%/M, 1/%/M), as 
for HypT-unreg and Sigm-unreg cases, input weights are automatically kept small 
when the network size increases, thus exploitin g the central part of bot h activation 
functions: consequently saturation is avoided ( Cancelliere et al I l2ni2|) . 

After clarification of the main issues affecting unregularised approaches, we 
now discuss the regularised one, derived according to eq. ®. 
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Fig. 3 Performance comparison in the Sigm (-reg and -unreg) cases; Abalone (Left); Iris 
(Right) 


Looking at HypT-reg and Sigm-reg cases in Figures [2] and |3l it appears that 
not only numerical instability (i.e. the error peak) is removed, but also that the 
penalty term provides control of overfitting, avoiding error growth and allowing 
optimal exploitation of the superior potential of larger architectures. The error 
curves feature now a monotonic decrease, becoming increasingly smoother. 

We obtained this result basing on the implications of eq. m that clearly sug¬ 
gests the role of the A parameter in preventing instability: our original idea is to 
address the issue of its determination with an ‘ad hoc’ tuning whenever growing 
error drift is experienced by the unregularised approach, as follows. 

We evaluate the validation error trend inside the critical region, looking for its 
minimum as a function of A; this allows to select a suitable value for this parameter 
for the subsequent experimentation with regularised pseudoinversion. 

We then made 100 random input weight choices and evaluated the mean test 
error (RMSE or misclassification rate) and the standard deviation S for any num¬ 
ber of hidden neurons. In Table[2]we list for HypT-reg and Sigm-reg the best ones 
of these mean values, called Optimal performance, together with the number of 
hidden neurons used to reach the associated performance (in parenthesis), and the 
corresponding standard deviation, as well as the selected value of A. 

Error values significantly better, basing on the Student’s t-test evaluation of 
statistical confidence intervals, are recorded in bold. We can see that there is 
roughly no “winner” between HypT-reg and Sigm-reg. 

Some other interesting considerations can be made noting that the regularised 
error plots in Figures [2] and [S] after an initial rapid decrease, achieve a nearly 
constant regime, with small variation vs. increasing numbers of hidden neurons. We 
can thus dehne a “near optimal” network size in both cases, as the one associated 
to a near optimal error level, i.e. an error signihcantly better than that from 
other methods. The assessment is again based on the Student’s t-test evaluation 
of statistical conhdence intervals. 

Table |3] compares the performance of this near optimal network with those 
obtained with the other methods, listing the best mean error (with the associ¬ 
ated number of hidden neurons in parenthesis) and the corresponding standard 
deviation. 

Thus, it appears that regularisation provides, except for Landsat dataset, 
the best performance not only in terms of lowest error values (see table [2]) but 
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Table 2 Optimal performance for regression (top) and classification (bottom) datasets. For 
Delta Ailerons, average errors and standard deviations have to be multiplied by 10“'^ 


HypT-reg 

Sigm-reg 

Abalone 

RMSE 

2.168 (192) 
2.150 (240) 

S 

0.003 

0.004 

A 

10-11 

10-12 

Cpu 

RMSE 

57.0 (78) 
57.3 (88) 

S 

2.1 

1.8 

A 

10-11 

10-10 

Delta Ail. 

RMSE 

1.639 (147) 
1.636 (244) 

S 

1.9e-3 

2.0e-3 

A 

10-13 

10-13 


Iris 



Diabetes 



Landsat 




Err.(%) 

S(%) 

A 

Err.(%) 

S(%) 

A 

Err.(%) 

S(%) 

A 

HypT-reg 

0.02 (66) 

0.2 

10-14 

20.2 (347) 

0.3 

10-11 

12.65 (486) 

0.2 

10-10 

Sigm-reg 

0.2 (19) 

1.1 

10-14 

19.4 (95) 

0.6 

10-12 

13.77 (396) 

0.3 

10-11 


Table 3 Method comparison for regression (top) and classification (bottom) tasks. For Delta 
Ailerons, average errors and standard deviations have to be multiplied by 10“"* 


HypT-reg (near opt.) 
Sigm-reg (near opt.) 

Abalone 

RMSE 

2.181 (32) 
2.181 (74) 

S 

0.011 

0.008 

Cpu 

RMSE 

57.55 (22) 
57.05 (45) 

S 

6.2 

3.2 

Delta Ail. 

RMSE 

1.642 (62) 
1.647 (99) 

S 

3.7e-3 

3.2e-3 

HypT-unreg 


2.187 (34) 

0.012 

59.37 (17) 

7.2 

1.649 (56) 

2.8e-4 

Sigm-unreg 


2.183 (32) 

0.013 

57.58 (16) 

8.9 

1.648 (45) 

6.5e-3 

ELM 


2.186 (33) 

0.015 

59.48 (19) 

9.1 

1.649 (75) 

8.8e-3 



Iris 


Diabetes 


Landsat 




Err.(%) 

S(%) 

Err.(%) 

S(%) 

Err.(%) 

S(%) 

HypT-reg (near opt.) 

0.56 (37) 

0.9 

21.0 (146) 

0.6 

12.65 (486) 

0.2 

Sigm-reg (near 

opt.) 

0.3 (16) 

1.1 

20.7 (45) 

0.8 

13.77 (396) 

0.3 

HypT-unreg 


1.02 (30) 

1.1 

21.2 (80) 

1.2 

12.94 (295) 

0.4 

Sigm-unreg 


1.00 (29) 

1.2 

20.8 (86) 

1.4 

12.78 (377) 

0.4 

ELM 


2.52 (23) 

1.9 

21.2 (69) 

1.2 

10.76 (390) 

0.4 


even with usage of smaller networks, limiting in this way computational cost and 
model complexity, an d therefore fulhlling the goals set by previous researches (e.g. 
( Yu and Pend . 120121) ). The smaller standard deviations almost always associated 
with the regularised methods also suggest a lower dependence from initial condi¬ 
tions. 

We also highlight that the use of small input weights and sigmoidal activation 
functions, which characterizes the Sigm-unreg case, allows to obtain error values 
lower with respect to the ELM case, so conhrming the effectiveness of this choice 
in order to contain saturation and overfitting issues. 

Landsat dataset constitutes an exception because best performance is reached 
using ELM. In this case, regularisation does not seem to be required, because 
overfitting and/or numerical instability do not take place, as it appears from un¬ 
regularised error curves. The behaviour of the threshold ratio, which remains al¬ 
ways larger than unity, is consistent with the lack of numerical instability and 
with our hypothesis of its relationship with error peaks. The lack of overhtting 
appears to be specific to the complexity of the dataset, having input vectors with 
size much larger than others, and therefore requiring a significantly larger number 
of parameters. 

In table |3] are listed the computational times (in seconds) recorded for all 
datasets for completing one training step, i.e. one random setting of input weights 
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Table 4 Comparison of training times (s.) for regression (top) and classification (bottom) 
tasks at fixed hidden layer size (100 neurons). 



Abalone 

Cpu 

Delta Ail. 

HypT-reg 

0.032 

0.035 

0.093 

Sigm-reg 

0.029 

0.040 

0.096 

HypT-unreg 

0.018 

0.018 

0.064 

Sigm-unreg 

0.0151 

0.025 

0.077 

ELM 

0.0123 

0.015 

0.060 


Iris 

Diabetes 

Landsat 

HypT-reg 

0.037 

0.0158 

16.65 

Sigm-reg 

0.032 

0.0155 

16.77 

HypT-unreg 

0.017 

0.0139 

15.94 

Sigm-unreg 

0.021 

0.0144 

15.78 

ELM 

0.011 

0.0116 

14.32 


and one output weights evaluation through pseudoinversion of the hidden layer 
output matrix; the number of hidden neurons has been fixed to 100 for the sake 
of comparison. 

We can see that, for each dataset, the times associated to each method do not 
differ significantly, because after having fixed the number of hidden neurons, the 
computational load necessary for the processing is comparable. 

The interested reader can find, for the com mon datatsets , a com parison among 
training times of ELM and backpropagation in ( Huang et al. I l2006l) . and can verify 
that ELM turns out to be two or three orders of magnitude faster. 


5 Conclusions 


We have considered the numerical instability ad overfitting problems for single 
hidden layer neural networks trained by pseudoinversion. We have shown how 
to use singular value analysis for the diagnosis of numerical instability, and how 
to solve this problem through the determination of a critical hidden layer region 
from which the regularisation technique benefits. This method also contributes to 
reduce the overfitting. 

Tests have been performed for both regression and classification tasks. For five 
out of six cases, the proposed regularisation is proven necessary and provides a 
significant performance improvement with respect to unregularised techniques; it 
also allows to built lean architectures which achieve near optimal performance 
with a reduced number of hidden neurons. 

Moreover the use of sigmoidal activation functions and “small” input weights 
(small because their values are linked to the hidden layer size), which characterizes 
the Sigm-unreg case, allows to obtain error values lower with respect to the ELM 
case, so confirming the effectiveness of this choice in order to contain saturation 
and overfitting issues. 

Comparing our results on the c ommon regression datasets with the alternative 
method proposed by Miche et al. ( Miche et al. IM), we note that our technique 
achieves RMSE values lower than those corresponding to their MSE values, with a 













An analysis of numerical issues in neural training by pseudoinversion 


11 


somewhat lower number of neurons. Besides, in our opinion, our method is simpler, 
in the sense that it uses a single step of regularisation rather than two in their 
method, and we also deal with classification tasks. 
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